
#include "gsl_ext_function.h"
#include "gsl_ext_multider.h"
#include "gsl_ext_symbolic_function.h"
#include "gsl_ext_testfunctions.h"

typedef struct
{
  double (*f)(const gsl_vector *x, void *params);
  void (*df)(const gsl_vector *x, void *params, gsl_vector *df);
  void (*d2f)(const gsl_vector *x, void *params, gsl_matrix *d2f);
  int n;
  gsl_ext_multimin_function_params *p;
}
func_alloc_data;

static int function_alloc(char *expr,
                          gsl_ext_func_eval_type fe,
                          gsl_ext_deriv_eval_type ge,
                          gsl_ext_deriv_eval_type he,
                          void *params,
                          func_alloc_data *result)
{
  result->p = malloc(sizeof(gsl_ext_multimin_function_params));

  result->p->f_eval_type = fe;
  result->p->df_eval_type = ge;
  result->p->d2f_eval_type = he;
  result->p->params = params;
  result->p->f_eval_counter = 0;
  result->p->df_eval_counter = 0;
  result->p->d2f_eval_counter = 0;

  if(fe == GSL_EXT_FUN_EVAL_SYM)
  {
    result->f = result->p->f = &gsl_ext_multimin_symbolic_function_f;
    result->p->f_evaluator = gsl_ext_multimin_symbolic_function_f_evaluator_alloc(expr);
    if(result->p->f_evaluator == NULL)
    {
      free(result->p);
      GSL_ERROR_VAL("invalid function expression", GSL_EBADFUNC, GSL_FAILURE);
    }
    result->n = ((gsl_ext_multimin_symbolic_function_f_evaluator *)result->p->f_evaluator)->n;

    if(ge == GSL_EXT_DERIV_EVAL_SYM)
    {
      result->df = result->p->df = &gsl_ext_multimin_symbolic_function_df;
      result->p->df_evaluator = gsl_ext_multimin_symbolic_function_df_evaluator_alloc(result->p->f_evaluator);

      if(he == GSL_EXT_DERIV_EVAL_SYM)
      {
        result->d2f = result->p->d2f = &gsl_ext_multimin_symbolic_function_d2f;
        result->p->d2f_evaluator = gsl_ext_multimin_symbolic_function_d2f_evaluator_alloc(result->p->df_evaluator);
      }
      else if(he == GSL_EXT_DERIV_EVAL_FD)
      {
        result->d2f = result->p->d2f = &gsl_ext_multider_forward_difference_d2f;
        result->p->d2f_evaluator = gsl_ext_multider_evaluator_alloc(result->n, result->p);
      }
    }
    else if(ge == GSL_EXT_DERIV_EVAL_FD)
    {
      result->df = result->p->df = &gsl_ext_multider_central_difference_df;
      result->p->df_evaluator = gsl_ext_multider_evaluator_alloc(result->n, result->p);

      result->d2f = result->p->d2f = &gsl_ext_multider_forward_difference_d2f;
      result->p->d2f_evaluator = result->p->df_evaluator;
    }
  }
  else
  {
    /* finite-difference gradient must be used 
     * if symbolic function expression is not available */
    if(ge == GSL_EXT_DERIV_EVAL_SYM)
      GSL_ERROR_VAL("Symbolic gradient requires symbolic function expression.", GSL_EINVAL, GSL_FAILURE);
    /* finite-difference hessian must be used 
     * if symbolic expression is not available */
    if(he == GSL_EXT_DERIV_EVAL_SYM)
      GSL_ERROR_VAL("Symbolic Hessian requires symbolic function expression.", GSL_EINVAL, GSL_FAILURE);
    const gsl_ext_testfunction_type *tf = gsl_ext_get_testfunction(expr);
    /* the given function name does not match 
     * any predefined test function */
    if(tf == NULL)
      GSL_ERROR_VAL("The given function name does not match any predefined test function.", GSL_EINVAL, GSL_FAILURE);
    /* no parameters supplied to a test function 
     * that requires parameters, e.g. dimensions */
    if(params == NULL && gsl_ext_testfunction_needs_params(tf))
      GSL_ERROR_VAL("The requested test function requires parameters.", GSL_EINVAL, GSL_FAILURE);

    if(params != NULL)
    {
      gsl_ext_testfunction_params *tfp = (gsl_ext_testfunction_params *)params;
      if(!gsl_ext_testfunction_check_params(tf, tfp))
        GSL_ERROR_VAL("Invalid test function parameters given.", GSL_EINVAL, GSL_FAILURE);
      if(tf->n == 0)
        result->n = tfp->n;
      else
        result->n = tf->n;
    }
    else
      result->n = tf->n;

    result->f = result->p->f = tf->f;
    
    if(ge != -1)
    {
      result->p->df_evaluator = gsl_ext_multider_evaluator_alloc(result->n, result->p);
      result->df = result->p->df = &gsl_ext_multider_central_difference_df;
    }
    if(he != -1)
    {
      result->p->d2f_evaluator = gsl_ext_multider_evaluator_alloc(result->n, result->p);
      result->d2f = result->p->d2f = &gsl_ext_multider_forward_difference_d2f;
    }
  }

  return GSL_SUCCESS;
}

gsl_multimin_function *gsl_ext_multimin_function_alloc(char *expr,
                                                       gsl_ext_func_eval_type fe,
                                                       void *params)
{
  func_alloc_data alloc_results;
  gsl_multimin_function *f;

  if(function_alloc(expr, fe, -1, -1, params, &alloc_results) != GSL_SUCCESS)
    GSL_ERROR_VAL("unable to initialize the user-supplied function", GSL_EBADFUNC, 0);

  f = malloc(sizeof(gsl_multimin_function));
  f->f = alloc_results.f;
  f->params = alloc_results.p;
  f->n = alloc_results.n;

  return f;
}

void gsl_ext_multimin_function_free(gsl_multimin_function *f)
{
  gsl_ext_multimin_function_params *p = (gsl_ext_multimin_function_params *)f->params;

  if(p->f_eval_type == GSL_EXT_FUN_EVAL_SYM)
    gsl_ext_multimin_symbolic_function_f_evaluator_free(p->f_evaluator);

  free(f->params);
  free(f);
}

gsl_multimin_function_fdf *gsl_ext_multimin_function_fdf_alloc(char *expr,
                                                               gsl_ext_func_eval_type fe,
                                                               gsl_ext_deriv_eval_type ge,
                                                               void *params)
{
  func_alloc_data alloc_results;
  gsl_multimin_function_fdf *f;

  if(function_alloc(expr, fe, ge, -1, params, &alloc_results) != GSL_SUCCESS)
    GSL_ERROR_VAL("error with initializing the user-supplied function", GSL_EBADFUNC, 0);

  f = malloc(sizeof(gsl_multimin_function_fdf));
  f->f = alloc_results.f;
  f->df = alloc_results.df;
  f->fdf = &gsl_ext_multimin_fdf;
  f->params = alloc_results.p;
  f->n = alloc_results.n;

  return f;
}

void gsl_ext_multimin_function_fdf_free(gsl_multimin_function_fdf *f)
{
  gsl_ext_multimin_function_params *p = (gsl_ext_multimin_function_params *)f->params;
  
  if(p->df_eval_type == GSL_EXT_DERIV_EVAL_SYM)
    gsl_ext_multimin_symbolic_function_df_evaluator_free(p->df_evaluator);
  else
    gsl_ext_multider_evaluator_free(p->df_evaluator);
  if(p->f_eval_type == GSL_EXT_FUN_EVAL_SYM)
    gsl_ext_multimin_symbolic_function_f_evaluator_free(p->f_evaluator);
  
  free(f->params);
  free(f);
}

gsl_ext_multimin_function_fd2f *
gsl_ext_multimin_function_fd2f_alloc(char *expr,
                                   gsl_ext_func_eval_type fe,
                                   gsl_ext_deriv_eval_type ge,
                                   gsl_ext_deriv_eval_type he,
                                   void *params)
{
  func_alloc_data alloc_results;
  gsl_ext_multimin_function_fd2f *f;

  if(function_alloc(expr, fe, ge, he, params, &alloc_results) != GSL_SUCCESS)
    GSL_ERROR_VAL("error with initializing the user-supplied function", GSL_EBADFUNC, 0);

  f = malloc(sizeof(gsl_ext_multimin_function_fd2f));
  f->f = alloc_results.f;
  f->df = alloc_results.df;
  /*if(ge == GSL_EXT_DERIV_EVAL_SYM)*/
    f->fdf = &gsl_ext_multimin_fdf;
  /*else
    f->fdf = &gsl_ext_multider_forward_difference_fdf;*/
  f->d2f = alloc_results.d2f;
  /*if(he == GSL_EXT_DERIV_EVAL_SYM || ge == GSL_EXT_DERIV_EVAL_SYM)*/
    f->fd2f = &gsl_ext_multimin_fd2f;
  /*else
    f->fd2f = &gsl_ext_multider_forward_difference_fd2f;*/
  f->params = alloc_results.p;
  f->n = alloc_results.n;

  return f;
}

void gsl_ext_multimin_function_fd2f_free(gsl_ext_multimin_function_fd2f *f)
{
  gsl_ext_multimin_function_params *p = (gsl_ext_multimin_function_params *)f->params;

  if(p->d2f_eval_type == GSL_EXT_DERIV_EVAL_SYM)
    gsl_ext_multimin_symbolic_function_d2f_evaluator_free(p->d2f_evaluator);
  else
    gsl_ext_multider_evaluator_free(p->d2f_evaluator);
  if(p->df_eval_type == GSL_EXT_DERIV_EVAL_SYM)
    gsl_ext_multimin_symbolic_function_df_evaluator_free(p->df_evaluator);
  else
  {
    if(p->df_evaluator != p->d2f_evaluator)
      gsl_ext_multider_evaluator_free(p->df_evaluator);
  }
  if(p->f_eval_type == GSL_EXT_FUN_EVAL_SYM)
    gsl_ext_multimin_symbolic_function_f_evaluator_free(p->f_evaluator);

  free(f->params);
  free(f);
}

void 
gsl_ext_multimin_fdf(const gsl_vector *x, void *params, double *f, gsl_vector *g)
{
  gsl_ext_multimin_function_params *p = 
    (gsl_ext_multimin_function_params *)params;

  *f = (*p->f)(x, params);
  (*p->df)(x, params, g);
}

void 
gsl_ext_multimin_fd2f(const gsl_vector *x, void *params, double *f, gsl_vector *g, gsl_matrix *H)
{
  gsl_ext_multimin_function_params *p = 
    (gsl_ext_multimin_function_params *)params;

  *f = (*p->f)(x, params);
  (*p->df)(x, params, g);
  (*p->d2f)(x, params, H);
}
